Analysis of the Platynereis dumerilii connectome with CATMAID and R

Gáspár Jékely
Centre for Organismal Studies, Heidelberg University

Resources



Part 2 - Packages, project management, data import, tidy data


  • some rules on files, folders
  • R projects
  • installing and loading packages
  • project templates
  • importing various types of datasets
  • principles of tidy data, tidying of messy datasets

Go to GitHub and fork my repository



  • Clone the repository: https://github.com/JekelyLab/ELTE_Catmaid_course_25
  • go to RStudio -> New Project -> Git -> repo’s address
  • select local dir and download the project with all directories to your computer
  • let’s check the project
  • open /analysis/scripts/Course_exercises.R

Rule nr. 1 – relative working directories


getwd()
  • Never use absolute paths in your scripts, because they hinder sharing: no one else will have exactly the same directory configuration as you.

  • Do not use setwd() to set your working dir

Use R projects (.Rproj)


  • Keep all files associated with a project together — input data, R scripts, analytical results, figures. This is such a wise and common practice that RStudio has built-in support for this via Rprojects.

  • If you create a new Rproject, your working dir will in general be where you save the new project

  • Whenever you refer to a file with a relative path it will look for it in your working dir.

  • in case you used setwd() in a script (you should not) you can find your Rproject directory with here()

here::here()

#test, but you should not do this
setwd("~")
getwd()

#get back to Rproj dir
setwd(here::here())
getwd()
# all good again

Use R projects (.Rproj)


  • list files in your working dir
list.files()
 [1] "analysis"                            "CITATION.cff"                       
 [3] "custom.scss"                         "ELTE_CATMAID_Presentation_files"    
 [5] "ELTE_CATMAID_Presentation.html"      "ELTE_CATMAID_Presentation.qmd"      
 [7] "ELTE_CATMAID_Presentation.rmarkdown" "LICENSE"                            
 [9] "manuscript"                          "new_paper_template.Rproj"           
[11] "README.md"                           "README.Rmd"                         
[13] "sessionInfo.txt"                     "versionInfo.txt"                    

Project management




  • Use folders relative to your main .Rproject file (e.g. My_next_report.Rproject)

  • Use a consistent directory structure to store code, data, text, figures, supplements, etc.

  • Can be ensured if you always use the same project template

  • Check the folder structure in the downloaded project template

Save and share your computer environment and packages

#save session info and Rstudio version info for reproducibility
sessionInfo()
R version 4.4.1 (2024-06-14)
Platform: x86_64-pc-linux-gnu
Running under: Ubuntu 24.04.1 LTS

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.12.0 
LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.12.0

locale:
 [1] LC_CTYPE=en_GB.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_GB.UTF-8        LC_COLLATE=en_GB.UTF-8    
 [5] LC_MONETARY=en_GB.UTF-8    LC_MESSAGES=en_GB.UTF-8   
 [7] LC_PAPER=en_GB.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_GB.UTF-8 LC_IDENTIFICATION=C       

time zone: Europe/London
tzcode source: system (glibc)

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

loaded via a namespace (and not attached):
 [1] compiler_4.4.1    fastmap_1.2.0     cli_3.6.3         tools_4.4.1      
 [5] htmltools_0.5.8.1 rstudioapi_0.17.1 yaml_2.3.10       rmarkdown_2.29   
 [9] knitr_1.49        jsonlite_1.8.9    xfun_0.50         digest_0.6.37    
[13] rlang_1.1.5       evaluate_1.0.3   
writeLines(capture.output(sessionInfo()), "sessionInfo.txt")

Installing packages

  • install the tidyverse package
install.packages("tidyverse")
  • then load the package
library(tidyverse)
library(png)

Source packages and functions from one file



  • open the R script analysis/scripts/Course_exercises.R with the example code

  • source several packages and functions, listed in one file consistently across scripts

source("analysis/scripts/packages_and_functions.R")

Import neurons from CATMAID and plot

Code
#read the skeleton ids first
neurons_skids <- skids_by_3annotations(
  "episphere", "Sensory neuron", "adult eye"
  )

#load neurons
neurons <- nlapply(
    read.neurons.catmaid(
      neurons_skids,
      pid = 11
    ),
    function(x) smooth_neuron(x, sigma = 6000)
  )

plot_background_anterior()

plot3d(
  neurons,
  WithConnectors = F, soma = T, lwd = 2,
  add = T, alpha = 0.4,
  col = "red"
  )


# Convert the rgl plot into an interactive htmlwidget
rglwidget()

Retrieve synaptic partners

#find synaptic partners
partners <- catmaid_query_connected(
  neurons_skids, minimum_synapses = 10,
  pid = 11
  )
partners
$outgoing
     skid partner syn.count num_nodes
37  25010   32142        23      3009
2    7501   37580        20      2747
36  23921   32142        20      3009
1    6743   37580        17      2747
3   10536   37580        17      2747
38  26163   32142        17      3009
20  19515   22159        14      2230
59  14162   34044        14      2066
63  58694   34044        14      2066
34  13459   32142        12      3009
21  22876   22159        11      2230
58  13459   34044        11      2066
14 151847   37580        10      2747
39  27064   32142        10      3009

$incoming
[1] skid      partner   syn.count num_nodes
<0 rows> (or 0-length row.names)

Parse IDs and load synaptic partners

#parse partner IDs
outgoing_partner_skids <- unique(
  partners$outgoing$partner
  )

#load partner neurons
outgoing_partners <- nlapply(
  read.neurons.catmaid(
    outgoing_partner_skids,
    pid = 11
  ),
  function(x) smooth_neuron(x, sigma = 6000)
)

Read neuron names and plot partners

#get name of neurons
catmaid_get_neuronnames(
  outgoing_partner_skids, pid = 11
  )
   32142    37580    22159    34044 
"IN1_pl" "IN1_pr" "IN1_al" "IN1_ar" 
#plot the partners to the same scene
plot <- plot3d(
  outgoing_partners,
  WithConnectors = F, soma = T, lwd = 2,
  add = T, alpha = 1,
  col = blues[8]
)

# Convert the rgl plot into an interactive htmlwidget
rglwidget()

Add text label to RGL scene

#add text label
texts3d(
  35000, 40000, 21000, 
  text = "PRC", 
  col = "black", cex = 2
  )

texts3d(
  35000, 62000, 21000, 
  text = "IN1", 
  col = "black", cex = 2
)

rglwidget()

Save RGL scene as png image

#save rgl view as png
rgl.snapshot("analysis/pictures/eye_with_IN1.png")

#close RGL viewer
close3d()